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ABSTRACT 


Closed -form analytic expressions for the time variations of 
instantaneous orbital parameters and of the topocentric range and range rate 
of a spacecraft moving in the gravitational field of an oblate large body are 
derived using a first-order Variation of Parameters technique. In addition, 
the closed-form analytic expressions for the partial derivatives of the topo- 
centric range and range rate are obtained, with respect to the coefficient of 
the second harmonic of the potential of the central body (J^). The results 

are applied to the motion of a point-mass spacecraft moving in the orbit 

-5 

around the equatorially elliptic, oblate Sun, with ^ 2.7 X 10 
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I. INTRODUCTION 


The results of the theory of motion of a particle moving in a close 
orbit around a large oblate body are well known in the literature. Although 
the equation of the trajectory cannot be obtained in the closed form, a multi- 
tude of solutions to this problem have been found in many different forms, 
with different orders of accuracy, by means of power series of the coeffi- 
cient of the second harmonic J^. 

The ultimate purpose of this study is twofold: first to determine the 
effects of the additional disturbing force generated by the oblateness of the 
central body on the observables — topocentric distance (range) and radial 
velocity (range rate) of a point -mass spacecraft-, and secondly, to furnish 
the analytic expressions for partial derivatives of the two observables 
needed for the trajectory improvement analysis. 

For close-to-Sun missions, assuming an oblate Sun having an oblateness 

_ 5 

coefficient = 2. 7 X 10 , due to the fact that the coefficient is very small, 

a first-order perturbation theory would yield sufficiently accurate results. 

II. POTENTIAL FUNCTION OF THE CENTRAL BODY AND 
THE EQUATION OF MOTION OF THE SPACECRAFT 

The potential of the attraction force exerted upon a particle by a large 
oblate body of mass M is given by (Refs. 1, 2) 

v = 

r 

where r is the magnitude of the position vector of the particle relative to 
the center of mass of the large body, are coefficients of zonal harmonics, 


® /F 

\ n 

■ ■ I >„(- 

n=2 

r) p „h> 
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Rq is the equatorial radius of the large central body, and GM is its 
gravitational constant; P^( 0 are Legendre polynomials of argument 

I* = sin 4> 

where 4> is the latitude of the particle measured from the equatorial plane 
of the central body. 

Considering only the equatorial ellipticity of the central body, we can 

- J 2 (t) P 2<^ (: 


write the potential in the form 


V = - 


GM 


where P^( £) is the Legendre polynomial of the second order 


P 2 U) = | (3; 2 - 1) 


The acceleration of the attraction force, having a potential given by Eq, (1), 
is 

a = = -grad V = -VV 

Hence, the acceleration can be written in the form 


dv _ GM- , - 

d t ‘ 3 r + a P 

r 


( 2 ) 


where ap is the additional, disturbing acceleration, caused by the equatorial 
ellipticity of the large body, 



■ GMJ 2 R 0 V 


P,U) 


(3) 


and r is the heliocentric position vector of the particle. 
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Let us introduce now an inertial frame of reference xyz, with the 
origin at the center of mass of the central body and the xy-plane coinciding 
with the equatorial plane of the central body. Then, 

Z* — — — 

sin <f> = t, = — , r = xi + y j + zk 


Performing the gradient operation indicated in Eq. (3), we find, for 
the disturbing part of the acceleration, 


a 


P 


! GMJ 2 




As it is customary in the theory of orbital perturbations, we shall 
assume the solution of Eq. (2) in the qua si'- elliptic form 


P 

1 + e cos 0 



where the semi-latus rectum p and the eccentricity e are functions of time. 
Such an ellipse, different in shape and dimensions at each point, is known 
and referred to in the literature as the oscillating ellipse, since the true 
trajectory of the particle and this ellipse have a common point and the same 
velocity at each instant during the motion of the particle (Refs, 3, 4). Using 
mathematical language, we can say that the true trajectory of the particle 
and the ellipse given by Eq. (4) have the tangency (or osculation) of the 
second order at every point during the motion. We must add here that the 
polar angle 0, also called the true anomaly of the particle, differs from the 
true anomaly of the unperturbed elliptic motion due to the action of the dis- 
turbing force. 

Denoting by 

r 2 

<=|gmj 2 -2. (5) 

p 


JPL Technical Memorandum 33-594 



we can write the expression for the disturbing acceleration in the form 


For the Sun, 




(5 sin 2 $ - 1 ) ^ 


2k sin <f> 


= 2. 7 x 10" 5 (Ref. 5). Using the values (Ref. 6) 

GM = 1. 327 X 10 U km 3 / s 2 
R q = 6. 980 x 10 5 km 


( 6 ) 


we find, with p given in astronomical units, 




0. 520 
4 
P 


X 10 


14 km/s 2 


(7) 


Since « «1, we can use the first-order Variation of Parameters method 

2 

within the accuracy of 0( € ). The parameters of the orbit (constants of 
integration of the Keplerian motion) which we shall use in this study are 

a = semimajor axis of the osculating ellipse 

e = its eccentricity 

£2 = longitude of the ascending node, i. e. , the angular distance 

between the positive direction of the x-axis and the line of nodes, 
in the equatorial plane of the central body 

i - angle of inclination of the orbital plane to the equatorial plane of 
the central body 

w = angular distance between the direction of the ascending node and 
the periapsis, in the orbital plane. This angle is called the argu- 
ment of periapsis. 

The first two elements, a and e, define the dimension and shape of 
the osculating ellipse. The three angles, £2, i, and go, are the three Euler 
angles which describe the position in space of the orbital plane and the 


4 
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orientation of the osculating ellipse in that plane. These angles also 
determine the motion of the plane under the action of the disturbing force. 

The equation of the Keplerian, undisturbed motion is a second-order 
vectorial differential equation, which, when integrated, contains six con- 
stants of integration. The choice of the sixth constant of integration, which 
is needed to introduce the time at which the particle is located at a certain 
position, is a matter of preference. In this study, we shall choose this 
sixth orbital parameter in the following manner. In Kepler's third law, 

n 2 a 3 = GM = (i (8) 

the quantity n represents the so-called mean motion, i. e. , the constant 
angular velocity of a particle which, moving in a circle, has the same 
orbital period as the particle moving in an ellipse. The quantity 

M(t) = n(t - T Q ) 

where T^ is the time of the periapsis passage, is called the mean anomaly of 
the particle at time t, located at point P(t) (Fig. 1). Let us assume that, at 
a certain initial epoch t^, we have an osculating ellipse, which we shall take 
as a reference orbit. According to the definition of the osculating orbit, if 
at time t = t^ the perturbations stopped acting on the particle, it would con- 
tinue to move in the initial reference ellipse and, at time t would be located 
at P*(t) (Fig. 1). The mean anomaly of the particle in that position is M*(t). 
We choose the sixth element t to be the difference (Ref. 3) 

t - M(t) - M*(t) (9) 

The orbital geometry and Euler's angles f2, i, and <*> are shown in 
Fig. 2. From the spherical triangle shown in Fig. 3, we find 

L, = sin 4> = sin i sin ( 0 + w) (10) 

Since our ultimate goal is to find the effect of the disturbing force on 
the topocentric range and range rate of the spacecraft, we have to find the 
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coordinates of the spacecraft in an Earth equatorial system XYZ, where 
the XY -plane is the equatorial plane of the Earth. The geometry is shown 
in Fig. 4. If a and 6 are the right ascension and declination of the north 
pole of the Sun, then, from Fig. 4, 


VN = a - 270 deg 
i N = 90 deg - 6 

Taking the x-axis in the direction of the node N, we see that, in order 
to bring the coordinate system XYZ into the position of the system xyz, we 
have to perform the following two rotations: 

(1) A positive rotation about the Z-axis by the angle YN; this rota- 
tion brings the X-axis into the position of the x-axis; 

(2) A positive rotation about the x-axis by the angle i^. This rota- 
tion brings the Z-axis into the position of the z-axis. 

Therefore, we can write the transformation equations 


and vice-versa, 


where 



X ' 


X 

Y 

“MINI 

y 

Z 


z 

r * 

X 


X ' 

y 

= M x N z 

Y 

z 


L z 


r* ° 

0 

- 

0 sin 6 cos 6 

X 

0 -cos 6 t 

sin 6 


( 11 ) 


( 12 ) 


(13) 
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and 


-sin a 



-cos a 
0 


cos a 0 

-sin a 0 

Q 1 


( 14 ) 


With respect to the mean equinox and Earth equator of 1950.0, the 
spherical coordinates of the north pole of the Sun are (Ref. 7) 


a = 286 S’ 0193 1 
6 = 6 3 f 7 7 1 8 ( 


(15) 


The three Euler angles i, £2, and u> of the orbital plane of the spacecraft, with 
respect to the equatorial plane of the central body (Sun), can be calculated 
from the corresponding angles i', £2', and w', given with respect to the equa- 
torial plane of the Earth, in the following manner. From the spherical 
triangle NNjN' (Fig. 4), we find 

cos i = sin 6 cos i' - cos 6 sin i ' sin (a - £2') 

sin i cos (u>' -w) = sin 6 sin i 1 + cos 6 cos i 1 sin (a - £2') 

sin i sin (u 1 - w) = -cos 6 cos (a - £2') 



sin i cos £2 = -cos 6 cos i' - sin 6 sin i' sin (a — £2 ' ) 


sin i sin £2 = -sin i 1 cos (a - £2') 


/ 


III. TWO INTEGRALS OF MOTION 


The differential equation of motion (2), with the disturbing accelera- 
tion ap given by Eq. (6), admits two integrals: the energy (vis -viva or 
Laplace's) integral and one scalar angular momentum integral. Indeed, 
writing the equation of motion in the form 
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W = - grad v 


and dot -multiplying both sides of this equation by vdt = dr, we find 


v . dv = -dr • grad V 


and, since dr • grad V = dV, the integration of the above equation yields 


2 v~ 

r= -v + vto) 


or 


= 2^-f £ p(|) 3 p 2 U) + H 


(17) 


where the constant H has the value 


tt 2 2 _ GM 

H = v. + 2V(0) = v A - 2 + 

0 0 r Q 


!<p(^) 3p 2<4o> 


Since for the elliptic osculating motion at t^ = 0, 


Vq = GM 


(n-A) 


the last equation for H transforms into 


H = 


2E 0 + !*p(^) P 2<V 


(18) 


where 


E 0 = 


GM 

2a(0) 


is the total energy of the particle at t^ = 0; a(0) is the value of the semi- 
major axis at that instant. Also, r^ = r(0) and = £(0). 
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To obtain the scalar angular momentum integral, we shall cross- 
multiply the differential equation (2) by vector r from the left. Thus, we 
obtain 


dh 

dt 


r x a p 


where h = r x v is the angular momentum vector. Substituting the value of 
a p from Eq. (6), we o,btain 


dh 

dt 


4 _ _ 

" 2 € ^(t) (r x k) 


(19) 


Vector r x k lies in the equatorial plane of the central body. Hence, 


dh __ _ 

= ( r x k ) • k = 0 


where h is the component of the vector h along the z-axis. Integrating the 
z ~ 

last equation, we have 


h = constant = h (0) 
z z 


From Fig. 2, we see that 


h = h * k = h cos i 
z 


and the scalar angular momentum integral becomes 


h cos i = h^ cos ip. 


i Q = 1(0), h Q = h(0) 


( 20 ) 
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IV. PERTURBATIONS OF THE OSCULATING 
ORBITAL PARAMETERS 


In order to obtain the time variations of the osculating orbital param- 
eters due to the effect of the disturbing force, we shall first utilize the 
already obtained integrals of motion given by Eqs. (17) and (20). Since, for 
the oscillating ellipse, 


v = GM 


( 1-0 


we obtain by subtraction 


1 1 


a(0) a 


2J 2 R 0 


p 2 (U P 2 (i 0 > 


Further, 


__1 J _ 

a(0) a 


a - a(0) 
a(0)a 


a - a <°L o(j 2 ) 

[a(0)] 2 V Z > 


and because of «.l, we can write 


Aa = 


-it 


' 3p ( a O P o) (?) P 2 U) _ t = 


where a^ = a(0), p n = p(0), and Aa = a - a(0). 


0 


( 21 ) 


( 22 ) 


Taking the true anomaly as the variable, we obtain from Eq. (4) 


( t ) 3 = 0 + -^) +3e (‘ +4-) coee 

Also, from 


cos 20 + -jp cos 36 


(23) 


P^f ?,) = y [3 sin^i sin^ (0+w) - 1 ] 


10 
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we find 


P 2 ^) ~ \ sin^ i - 2) - 3 sin^ i cos (20 + 2 oj)J 


(24) 


If we multiply Eqs. (23) and (24) and express everything in terms of multiple 
angles, we obtain, instead of Eq. (22), the following expression: 


Aa = J^K(a) 


(2-3 sin^ i) j (4 + e ^) cos 0 + 2e cos 20 + cos 3 0 


+ 


. 2 . 
sin i 


cos ( 6 - 2u>) + 3e(4 + e^) cos ( 0 + 2oo) 


+ 4(2 + 3e ) cos (20 + 2 uj) + 3e(4 + e ) cos (30 + 2 u>) 


+ 6e^ cos (40 + 2u>) + e^ cos (5 0 + 2u>) J 


1 G 

-e 


o 


(25) 


where 0^ = 0(0). The quantities e, i, and on the right-hand side of Eq. (25) 
should be written as e^, i^, and but for the sake of simplicity, we shall 

drop the subscript zero in further writing. The constant K(a) is given by 


v( ) _ 3 (V) 
K(a) - 8 T~ 


( 26 ) 


where the subscript zero has been dropped again. 

The second integral of motion, given by Eq. (20), yields 


di _ cot 
dt h 


cot i dh cot i /— dh \ 
h dt ■ h 2 \ h * dry 


where h = |h|. Substituting in this equation dh/dt given by Eq. 19), 
we find 
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di 

dt 


- 2 € ; 


cot 


m 


( r x k ) • h 


For the instantaneous osculating orbit, 


2 d6 
r dt 


= h 


( 27 ) 


so that, taking the true anomaly 0 as the independent variable, we obtain 

H = - 3 h(^-) cosi(|)sin ( e +U )|i- (28) 

We shall now substitute p/r from Eq. (4) and 

j y 

= sin i cos (0+w) 
d0 

which we obtain by differentiating Eq. (10) with respect to 0 . (For this 
differentiation, we keep i and w constant. ) After the substitution, Eq. (28) 
becomes 


where 


di 


= -6J 2 'K(i)(l + e cos0) sin (20 + 2u>) 


i/M 

8 \p / 


K (i) = q — J sin 2i 


(29) 


(30) 


Rewriting Eq. (29) i-n the form 


di 

d0 


- 3J 2 K(i) [ e sin ( 0 + 2u) + 2 sin (20 + 2co) + e sin (30 + 2 U )] 


(31) 


12 
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and integrating, keeping e, i, and u> fixed, we obtain 


Ai = J K(i) [3e cos ( 0 + 2«) + 3 cos (20 + 2w) + e cos (30 + 2w)l (32) 

L Je 0 

where 

Ai = i - i(0) 

To find the time variation of the angular distance of the ascending node, ft, 
we use the relationship (Ref. 3) 


flft _ tan ( 0 + u>) di 
d0 sin i d0 

Substituting di/d0 from Eq. (29), we obtain 

-4§- = - 6 J_ fl - cos (20 + 2w)"| ( 1 + e cos 0 ) 

d0 2 sin l L 'J 

After the integration, we find 

Aft = - J^K^ft) 6 ( © 4- e sin 0 ) - 3e sin ( 6 + 2w) 


where 


and 

Aft = ft - ft(0) 

As previously, the parameters e, i, and oo have been kept constant 
during this integration. 


0 


-3 sin (2 0+ 2w) - e sin (3 0 + 2 to) 

- e, 


(33) 


K(n) = ?nr) cos 1 


(34) 
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From Ref. 3 , we further find 


d oj . d£2 
dt + dt" cos>1 



where 6/dt is the well known variational operator. Taking the angle 0 
as the independent variable, the last equation can be written in the 
form 


dcj 

de 


df2 

de 


cos 


66 

d0 


( 35 ) 


Applying the variational operator to the polar equation of the osculating 
ellipse ( 4 ) and remembering that Sr = 0 * and that for every orbital parameter 


6 _ d 
dt ~ dt 


we find, taking the angle 6 as the independent variable, 

■ q 66 de Q 2p , . di 

e sin 6 — = -r— cos 6 tan l -3-3- 

d6 d0 r d6 

Also, from the definition of the semi-latus rectum, 



we find, by differentiation, 

de P da p , . di 

d6 ~ _ 2 dG ae an 1 de 
2a -e 


( 36 ) 


( 37 ) 


( 38 ) 


*The Brown-Milankovic variational operator 6/dt is defined in the two 
basic equations of the Calculus of Perturbations 


| 6 r = 0 
( 6 v = apdt 


14 
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Considering a, e, and i constant, we obtain, by integration, 


Ae = — - tan i 

ae V 2a 


i Ai^ 


( 39 ) 


Combining Eqs. (36) and (38), we find 


„ D 60 _ P cose da 2e + (1 + e ) cos 6 . di 

e sin 6 He z de i tan 1 He 

2a e 


(40) 


From Eq, (25), we obtain, by differentiation, 


= - J^K(a) £(2 - 3 sin 2 i) j (4 + e 2 ) sin 0 + 4e sin 26 + e 2 sin 30 j 
+ * |e 8 ^ ® " 2ui) + 3e(4 + e 2 ) sin ( 0 + 2u>) 


+ 8(2 + 3e 2 ) sin (20 + 2w) 


+ 9e(4 + e 2 ) sin (30 + 2c o) + 24e 2 sin (40 + 2w) 


+ 5e 2 sin (50' + 2co) | J 


(41) 


Subsituting expressions (41) and (31) into Eq. (40) and dividing the result by 
e sin0, we obtain 


60 _ 
d0 


J 2 K(w) 


(2-3 sin 


in 2 i ) ^ 


1 + - ^ g e - - cos 0 + cos 2 0 + cos 30 


) 


2 

+ — ^ — — |e cos (0 - 2co) - — — J e - cos (0 + 2w) 


+ 24 cos (20 + 2co) + 28 + cos (36 + 2w) 


4- 24 cos (40 + 2w) + 5e cos (56 + 2w) 


)] 


(42) 
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with 


KM = | 



(43) 


An integration of Eq. (35) with the angle i considered as constant yields 


Aw + Af2 cos i 



d 0 


or, after the performed integration on the right-hand side of the last 
equation. 


Aw + A£2 cos i = J^K(w) 


iin 2 i) ^0 


A + l 

(2-3 sin*" i) [0 + ^ — sin 0 + sin 20 + 


fz sin 3e ) 


. 2 . 
sin l 

24 


3e sin ( 0 - 2w) - — - 2 *^ - s in ( q 4 . 2w) 


+ 36 sin (26 + 2w) + 28 * lle sin (36 + 2w) 


+ 18 sin (46 + 2w) + 3e sin (50 + 2w) 


ie 

0 


o 


(44) 


Incidentally, it is interesting to note at this stage that the expressions 
for Aa, Ae, Ai, Af2, and Aw, given by Eqs. (25), (39), (32), (33), and (44), 
respectively, well agree with the known results. We find, namely, that Aa, 
Ae, and Ai are purely periodic, while Af 2 and Aw contain secular terms. 
Denoting by 



d 6 


16 


JPL Technical Memorandum 33-594 



the average value of an orbital element | over the period, we find that 


Hence, for 


[Aa] = 0, 


[Ae] = 0, 

C Ai] = 0, 

/ R \ ^ 

[A£I] = -f cos i, 

C Aw] = f J 2 ('t) (4 ' 5 Sin2 4) 


\ 


( 45 ) 


i 



there is no secular advancement of the periapsis line. 

According to what was said previously, the sixth orbital parameter, t , 
is the difference between the perturbed and unperturbed mean anomaly at any 
instant t. Since (Ref. 3) 


dM(t) _ , 6M(t) 

dt dt 


where 6/dt is the already used variational operator, and 


dM*(t) 
dt n 0 


we find by subtraction (see Eq. 9) that 


d_T 

dt 


n - n 0 + 


6M(t) 

dt 
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or 


d_T 

dt 


(46) 


“ = An + 


6M(t) 

dt 


Equation (46) can be written in a somewhat different form. From Kepler's 
third law (8), we obtain, by differentiation, 


An 


3n 

2a 


Aa 


so that Eq. (46) takes the form 


dT 3n A 6M(t) 

dt = - 2T Aa + ST 


(47) 


Kepler's equation 


4 - e sin 4 = M(t) (48) 

and the two relations which we derive from the Keplerian motion 

r cos 0 = a (cos 4 ~ e ) 


r sin 0 = b sin 4 


where 41 is the so-called eccentric anomaly angle, and 

b = a yjl - e ^ 


(49) 


is the semiminor axis of the osculating ellipse, enable us to obtain 6M(t)/dt 
in terms of 60/dt and de/dt. The expression which we find by using the 
operator 6/dt is 


6M(t) 

dt 



sin 0 


] 


18 
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or, taking again the true anomaly 0 as the independent variable, 


6M(t) 

d0 


2 r 

= (1 - e ) || - (2 + e cos 0) sin 0 e 


d0 


] 


( 50 ) 


Combining Eqs. (36), (38), and (50), we find 


6M(t) 

d0 


abe sin 


in 0 ( 


Pcos 0 - 2er da „ , 

- - p cos 0 tan i 


2a 


• ii\ 
1 d0 J 


(51) 


and, substituting (da)/(d0) and (di) / (d 0), respectively, from Eqs. (41) and 
(32), we obtain 


or p u 

6M(t) _ 2^0 

d0 4bpe 


(2-3 sin 2 i) 


M 1 '^) 


cos 0 - e cos 20 cos 30 

4 


. 2 . , ? 

+ — -|-e cos (0 - 2c o) + (4 + 17e ) cos (0 + 2 oj) 


+ 24e cos (20 + 2co) - (28 - 13e ) cos (30 + 2oo) 


- 24e cos (40 + 2u>) - 5e cos (50 + 2u>) 


(52) 


Although it is relatively simple to integrate Eq. (5 2), we shall do it by 
partially utilizing the already derived expression for Aco + Aft cos i in 
Eq. (44). Namely, subtracting expression (42) multiplied by V 1 - e^ from 
Eq. (52), we find 

- Vl - e 2 || = J 2 Vl - e 2 K(w)[2(2 - 3 sin 2 i)(l + e cos 0) 


+ 3 sin 2 i |e cos (0 + 2 cj) + .2 cos (20 + 2co) 

+ e cos (30 + 2w) J J (53) 
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Integrating the last equation, we obtain 


de + (Aw + AQ COS i)yjl - e 2 


J 2 K(w) ^2(2 - 3 sin 2 i)(0 + e sin 0) 

+ sin 2 i |3e sin (0 + 2w) + 3 sin (20 + 2w) 


+ e sin (30 + 2w) J 


-i0 

e, 


;0 


(54) 


Equation (47) in the form (h = nab) 


ix = .^_ (r2Aa) + aM!ti 


d0 


2a b 


yields, after the integration, 


At = - 


2a 2 b 


J (r 2 <Sa) de + J (^)de 


(55) 


where the last term on the right-hand side is given by Eq. (54), 
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From Eqs. (22) and (24), 


Aa = ^aR^ 2 [ 2(2 - 3 sin 2 i)(l + e cos 0) 


+ sin 2 i|3e cos (0 + 2u>) + 6 cos (20 + 2w) + 3e cos (30 + 2w) } ] 


+ 2J- 


(aR \ 2 

* p ( r ) 

3 *Z K *V 


(56) 


Hence, by integration. 


2a 2 b 


l 


(r Aa)d0 = 3J. 


naR ‘ 


o' 


At 


+~~ yj\ - e 2 K(w) 


2(2 - 3 sin 2 i)(0 + e sin 0) 


+ sin 2 i |3e sin (0 + 2w) + 3 sin (20 + 2w) 
+ e sin (30 + 2u>) 

u 0 



(57) 


Finally, substituting the results obtained in Eqs. (54) and (57) into Eq. (55), 
we find 
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At + (A 00 + Aft cos i) Vl - = - J^AAt 4- V 1 - e 2 K( W ) 


X (2-3 sin i)( 0 4- e sin 0) 


+ sin^ i 1 3 e sin (0 + 2u>) 4- 3 sin (20 4- 2u>) 


4- e sin (30 + 2u>) 


where 


3 3 P 2^0^ 


To find the change of t over one period, we write, from Eq. (58), 

[ At ] + ([ Aoj] + [A£2] cos i) yjl - e 2 = - J^A + J ^{2 - 3 sin 2 i) ^ 1 - e 2 K(oo) 


or, using the last two of Eqs. (45), 


[Ar]= -J 2 A = 


- 3J 2"^- P 2^0» 


From 


2 2 

3 sin i sin (0 q + w ) 
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it follows that, 


within the accuracy of O 




0 


if 0 


0 


is so chosen that 


provided that 


sin (0 + oj ) = ±- 


1 


\/3 sin i 


i > arc sin j = 35°3 
\\f3 / 

Equations (25), (32), (33), (39), (44), and (58) give the time variations 
of all six orbital parameters. 

V. SUMMARY OF FORMULAE FOR Aa, Ae, Ai, 

Af2, Agj, AND At 


A. Auxiliary Quantities 

We shall first compute the true anomaly angle 0 for the time t using 
the equation of the unperturbed motion. In other words, assuming no change 
in orbital parameters, we compute the eccentric anomaly angle 4 a t the 
time t from Eq. (48), 


4 - e sin 4 - n(t - t^) = nAt 


and the angle 0 from 


tan 0 


COS 4^ 


sin 4* 

- e 
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Next, we shall introduce the quantities 




24 



Having found all the necessary auxiliary quantities, we shall proceed 
with the calculation of the orbital parameter time variations. 

B. Time Variations of Orbital Parameters 

Substituting all the auxiliary quantities into Eqs. (25), (32), (33), (39), 
(44), and (58), we obtain 



Group C3] 
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VI. PERTURBATIONS OF THE HELIOCENTRIC DISTANCE AND 
THE POLAR ANGLE OF THE SPACECRAFT 


From Ref. 3, it follows that the disturbing effects on the heliocentric 
distance and the polar angle 0 (the true anomaly) and their derivatives 
(radial and angular velocities, r and 0 , respectively) are 
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In the first two equations of Group [4j, Ar and A0 are the differences 
between the heliocentric position and the polar angle of the particle in its 
perturbed trajectory at time t and the corresponding quantities related to 
an imaginary particle moving in the initial undisturbed elliptic orbit at the 
same time t. The same explanation holds for Ar and A0, given by the 
second two equations of Group [4]. 

VII. PARTIAL DERIVATIVES OF ORBITAL PARAMETERS 

WITH RESPECT TO J 2 

Since the orbital parameters a^, e^, i^, and M*(t) of the 

initial osculating reference orbit do not depend on J^, the partial derivative 
of any orbital parameter £ with respect to is 

di _ 9(A£) 

8J 2 - 8J 2 

Hence, we can write immediately 



JPL Technical Memorandum 33-594 


2.7 



= (DOJ) = -K(n)(6A0 + 6eS 1 - 3eL 1 - 3L 2 -eL 3 ) | 

2 


d i^ = (DWJ) = - (DOJ) cos i - K(u) 
9J 2 




I,1A9- + 4 + 3e S 
el 


2 S 2 + 12 S s) + 24 \ 


3e(L 0 + L 5 ) 


12- 21e 2 T , 28 + lie 2 T 

1 J— * -I + 0 

e 1 e 3 


+ 18(3L 2 + L 4 ) 


j T) = (DTJ) = - [(DWJ) + (DOJ) cos i]Vl - e 2 

- AAt + K(w)Vl " e 2 I (A0 + eL ) 


+ ^I 2 (3eL 1 + 3L 2 + eL 3 ) 


] 


Group C5] 


VIII. CALCULATION OF ORBITAL PARAMETERS i, ft, w FROM 
EQUATORIAL ORBITAL PARAMETERS i', ft', AND w' 


Most of the time in practice, the position of the orbital plane of a 
poi.nt-mass spacecraft, defined by the three Euler angles i, ft, and u> , is 
given relative to one of the Earth reference planes: the equatorial or ecliptic 
plane for a certain fixed epoch, say 1950. 0.- The orbital elements i, ft, and 
u> which we have used so far are relative to the equatorial plane of the central 
body, in particular, the equatorial plane of the Sun. To transfer from one 
system of Euler angles, i* , ft', w' (Earth equatorial), which are normally 
given, to the other system of Euler angles, i, ft, and u>, which we need 
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for the calculation of perturbations of the position of the orbital plane in 
space, we have to use Eqs. (16). With a and 6 given in Eqs. (15), the 
expressions from which we calculate i, f2, and w are 



Group [6J 


Let us denote by P the topocentric postion vector of the spacecraft 
relative to the observer's position at point P (Fig. 5), which is uniquely 
determined by the geocentric position vector |. If r and r^, are the helio- 
centric position vectors of the spacecraft and the Earth, respectively, from 
the diagram shown in Fig. 5, we can write 

P = r - r E - | (61) 

Differentiating Eq. (6l) with respect to time, we find 

p = r - r E - | (62) 

The disturbing effect of the equatorial ellipticity of the central body (Sun) 
affects not only the motion of the spacecraft but also the motion of the Earth. 
The position of the observer, however, defined by the geocentric position 
vector £, remains unchanged, i. e. , unaffected by the disturbing force. 

Denoting by r* and r^, the undisturbed position vectors of the space- 
craft and the Earth, respectively, in other words the heliocentric position 
vectors of points in which the spacecraft and the Earth would be at time t if 
the perturbing effect were turned off at the time t^, then, also 
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p* - r* - r£, - £ 


By subtraction of the last equation from Eq. (6l), we obtain 


Ap = p - p = ( r - r* 


- (r E- r E> 


= Ar - Ar, 


( 63 ) 


and, similarly, 


Ap = p - p* = ( r - r*) - ( r £ - r* ) 


= Ar - Ar 

J±< 

= Av - Av^ (64) 

iii 


where v and v.^, are, respectively, the velocities of the spacecraft and the 
Earth. 

Since 

p • Ap = p A p 

we find, dot -multiplying Eqs. (6 1 ) and (63), that 

• A p = "p [ r Ar + r E Ar E " ( r E + ¥) ' - r • Ar £ j 

or, because £ /p « 1, 

Ap = ^(rAr + r £ Ar E - r E • Ar - r • Ar E ) (65) 

where, within the same order of accuracy, 

P = £ r 2 + r|, + 2( r • r £ ) -J (66) 

From 

• 9 = P P 
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we obtain, by differentiation. 


Ap = 



A p + p • Ap-pAp) 


(67) 


and, substituting expressions (63) and (64) into the last equation, we find, 
with the same accuracy as before, 

Ap = ]5- ( v - V E ) ' ( Ar " Ar E ) + ( r - r E ) * ( Av - Av e ) - P A P j (68) 

where Ap is given by Eq. (65), 


Let X, Y, Z be the heliocentric rectangular coordinates of the space- 
craft in the Earth equatorial space-fixed system of reference axes, and X , 

E 

Y Z the heliocentric equatorial coordinates of the Earth in the same 
system. In the xyz-reference frame, with the equatorial plane of the central 
body (Sun) as the fundamental plane, the coordinates of the spacecraft and the 
Earth and their velocity components are, respectively, (Fig. 4 and Eq. 12), 
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The orbital plane parameters of the Earth, i^, S2 , and w , relative 

E ill E 

to the equatorial plane of the central body, can be easily determined in a 
manner similar to that used in determining the orbital plane elements of the 
spacecraft, using the equations of Group [6]. For the Earth, we should set 
£2^ = 0, and i^ = * E , where * E is the obliquity of the ecliptic. In reality, 
however, the orbital plane of the Earth (using astronomers’ jargon, we can 
call it true ecliptic, or instantaneous ecliptic plane) moves as a result of 
perturbations, and its position may be determined relative to the fixed 
Earth equatorial plane for the epoch 1950. 0. From Ref. 7, we find, for 
ecliptic elements of the orbital plane of the Earth, 

M 9 

i E = 0?0 13076t - 0 ? 000009 t 

= 174f 40956 - Of 24l66t + 0f00006 t 2 } (69) 

= 287f67097 + Of 56494t + 0f00009 t 2 t 

where the time is expressed in so-called Julian centuries, i.e., units of 
36, 525 days. Similarly, the obliquity of the ecliptic is (Ref. 7) 

i E = f E (t) = 23 ° 4457888616 - Of 0130141669 t 

- 0 f 09445 X 10“ 5 t 2 + Of 05000 X 10" 5 t 3 (70) 

The geometry of the orbital plane of the Earth is shown in Fig. 6. 

The orbital elements relative to the equatorial plane of the Earth for the 
epoch 1950. 0 are obtained from the spherical triangle YN^N^. Equations 
corresponding to Eqs. (16) are 
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cos i' 
E 


. . 1 i 

sin 1 ^ cos ft^ 

iij 


. * I . .M . _tl 

cos cos - sm sin cos ft^ 


cos i'' sin + sin i" cos e_ cos ft" 
E E E E E 


. . ' • ^ i . .11 . — n 

sin ij, sin ft^ = sin i sin ft^ 


sin i^ cos-(u)g - oo^) = sin i^ cos + cos sin e cos ft^l 


. . t . / I 1 I X 

sin i E sin (u> E - o)^) 


sin sin ^ 

E E 


( 71 ) 


Thus, we obtain 



Group [81 


Taking now the line of nodes as the principal reference axis in the 
orbital plane of the spacecraft (Fig. 2), we see that the coordinates of the 
spacecraft in the plane are (r cos u, r sin u, 0), where 


u = G + <*> 


(72) 


The rectangular equatorial coordinates x, y, z can be derived from the 
coordinates in the orbital plane by two rotations: first, a positive rotation 
about the line of nodes by the inclination angle i; second, a negative rotation 
about the polar axis of the central body by the angle ft. In other words. 
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X 

y 

X 




r cos u 
r sin u 
0 



y 


z 




r 

rii 

0 


\ 


} (73) 


/ 


For the Earth, we have two similar equations: 




- 

cos 

\ 

E 

^E 

■ Ml Nl 

E E 

r E sinu E 


.‘ E . 


0 



(74) 


X E 

m rpi rp 

r E 


*E 

■ M 2 M x m. 

r E U E 


*E 

J 


0 

1 


Since the equations for the Earth are obtained just by adding the subscript E, 
we shall write only the equations for the spacecraft and keep in mind that 
the same equations exist for the Earth. In expressions (73) and (74), 



are rotation matrices 
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H 


cos 12 
sin 12 
0 


- sin 12 
cos 12 
0 


0 

0 

1 


h: 


i 

0 

0 


cos 1 


sin i 


0 

-sin i 
cos i 


( 75 ) 


[■I 


cos u 
sin u 
0 


- sin u 
cos u 
0 


0 

0 

1 


/ 


Three rotation matrices can be written for the Earth by adding the 
subscript E to all quantities involved. 

Let us now write the expressions for the disturbing effects in range 
and range rate of the spacecraft, Ap and Ap, given by Eqs. (65) and (68), 
respectively, in the following forms: 
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We have, however, 


_9_ 

an 


_a_ 

di 


_d_ 

9u 


L 0 ] z L ^ J z 
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Au = A6 



X. PARTIAL DERIVATIVES OF RANGE AND RANGE RATE 


The partial derivatives of the topocentric range and range rate can be 
derived directly from the equations for the perturbational changes Ap and 
A p by simply replacing the operator A with in the equations of 

Groups [9] and [10]. Hence, from the equations of Group [10], we obtain 
the partial derivatives of the rectangular coordinates x, y, z and velocity 
components x, y, z, with respect to J^. Denoting by (Refs. 8 and 9) 


0 X 
aj 2 

= (DXJ) 

9x e 
9J 2 ' 

(DXEJ) 

ay 

aj_ 

= (DYJ) 

9 *e 

aj 0 

(DYEJ) 

2 


2 


dz 

dj 

- (DZ J) 

9z e _ 
aj 

(DZEJ) 


n 

•X ‘ 

(DDXJ) 

d E 
aj. 

(DDXEJ) 

Z 


2 


n 

•HJ 

(DDYJ) 

ay E 

aj 0 

(DDYEJ) 

Z 


2 


0Z 

aJ 2 = 

(DDZJ) 

9z e 
9J 2 “ 

(DDZEJ) 


\ 
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where 


8r 

(DRJ) = -qT 2 


8u 

(DU J) = 


, _ 8_r_ 

(DDR J) - qj 


(DDUJ) = gj 


2 



( 80 ) 


for the spacecraft, and 


8r E 

(DREJ) = 


8u E 

(EUEJ) = 


9i> E 

(DDREJ) = 


(DDUEJ) 


8u 
8 J. 



( 81 ) 


for the Earth. 
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Since r* and 6* do not depend on J^, we find from the equations of 
Group Q4] and Eq. (7 2) that 


(DR J) = - (DAJ) - a cos 0 (DEJ) + sin 0 (DTJ) 

a D 


(DUJ) = J + 0 sin 6 (DEJ) + (DT J) + (DWJ) 

2 

(DDR J) =. - sin 6 (DAJ) + nb(j-j sin 6 (DEJ) 

+ nae(0 cos 0 (DTJ) 

3 

(DDUJ) = -^y(DAJ) +^— (^cos 0 - e)(DEJ) 

2r ^ 7 


br 


2ne(0 sin 6 (DTJ) 


Group Cl 2] 


and similar expressions for the partial derivatives of r„, , u„, and u 

iL L E E 

with respect to J^. 

Finally, we can write 

\ 

(DXJ) - (DXEJ)l 
(DYJ) - (DYEJ) 

(DZJ) - (DZEJ) 


dP 

9J„ 




x T 


y -Yi 


z - z, 


9P 

9J 2 


/ 




y-y T 



'(DXJ) - (DXEJ) " 
(DYJ) - (DYEJ) 
(DZJ) - (DZEJ) 


(82) 


+ [ x - X E * - y E 2 



(DDXJ) - (DDXEJ) 
(DDYJ) - (DDYEJ) 
(DDZJ) - (DDZEJ) 



/ 


The computational sequence in the form of a flow chart is shown in Fig. 7. 
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It should be mentioned, in conclusion, that a much higher accuracy is 
obtained if, instead of taking initial orbital parameters a^, e^, i^, £2^, and 
cOq in the equations of Group [3D and thereafter, we repeat the calculations of 
Aa, Ae, Ai, A£2, Aw,, and At with the mean values of elements 


a 

m 



+ a) 



a 

a Q = a(0) 



e Q = e(0) 


G 

h 


h Q = h(0) 


e 

m 


2 (e 0 


+ e) 


4 - 

e o + — 


etc. 


NOMENCLATURE 

semimajor axis of the osculating ellipse at time t 

semimajor axis of the initial osculating reference ellipse 
at time t^ 

acceleration vector of the particle 
disturbing acceleration vector 
constant of integration 

semiminor axis of the osculating ellipse 

difference between cosines of angle N 0 at time t and 
time tp, for N = 1, 2, 3 

eccentricity of the osculating ellipse 

eccentricity of the initial osculating reference ellipse 

total energy of the particle at time t^ 

universal constant of gravitation 

magnitude of the angular momentum vector 

magnitude of the angular momentum vector at time t^ 

component of the angular momentum vector along the z-axis 
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v°> 


h 

H 

i 


*0 = i(0) 


E 


N 


I I* *2 


3 

J. 


k 

K, 


K. 


N 


K(a), K(i),1 
K(Q), K(w) j 


component of the angular momentum vector along the 
z-axis at time t^ 

angular momentum vector 

constant of integration 

angle of inclination of the orbital plane of the spacecraft 
to the equatorial plane of the central body- 

angle of inclination of the orbital plane of the spacecraft to 
the equatorial plane of the central body at time t^ 

angle of inclination of the orbital plane of the spacecraft to 
the equatorial plane of the Earth for 1950. 0 

angle of inclination of the orbital plane of the Earth to the 
equatorial plane of the central body 

angle of inclination of the orbital plane of the Earth to its 
equatorial plane for 1950. 0 

angle of inclination of the orbital plane of the Earth to the 
ecliptic plane for 1950.0 

angle of inclination of the equatorial plane of the central 
body to the equatorial plane of the Earth for 1950. 0 

unit vector along the x-axis 

constants depending on the angle of inclination i 
unit vector along the y-axis 

coefficient of the second harmonic of the potential of the 
central body 

unit vector along the z-axis 

difference between the cosines of angle 0 - 2u at time t 
and time t^ 

difference between the cosines of angle N0 + 2wat time t 
and time t^ for N = 1, 2, • • • , 5 

coefficients of expressions for the time variations of 
osculating orbital parameters 
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M 

M(t) 


M*(t) 


n 

n 0 = n (°) 
N 

P 

P 0 = P(0) 

P n U) 

P 2<U 

r 


r 


0 


r 


E 




r 



r 


r 


0 


difference between the sines of angle 0 - 2oo at time t and 
time tp 

difference between the sines of angle N0 + 2u> at time t and 
time t^ for N = 1, 2, • • • , 5 

mass of the central body 

mean anomaly angle of the particle in its true orbit at 
time t 

mean anomaly angle of an imaginary particle moving in the 
initial osculating reference orbit at time t 

mean motion of the particle in its orbit 

mean motion of the particle at time t^ 

subscript 

semi-latus rectum of the osculating ellipse at the time t 
semi-latus rectum of the initial osculating reference ellipse 
Legendre polynomial of order n 
Legendre polynomial of order two 

magnitude of the heliocentric position vector of the 
spacecraft 

magnitude of the heliocentric position vector of the space- 
craft at time t^ 

magnitude of the heliocentric position vector of the Earth 

magnitude of the heliocentric position vector of an 
imaginary particle moving in the initial osculating reference 
ellipse at time t 

radial velocity of the spacecraft 

radial velocity of the Earth 

heliocentric position vector of the spacecraft at time t 
heliocentric position vector of the spacecraft at time t^ 
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heliocentric position vector of an imaginary particle moving 
in the initial osculating reference ellipse at time t 

heliocentric position vector of the Earth at time t 

heliocentric position vector of an imaginary Earth moving 
in an initial osculating reference ellipse without 
perturbations 

velocity vector of the spacecraft 
velocity vector of the Earth 

velocity vector of an imaginary particle moving in the initial 
osculating ellipse without perturbations 

velocity vector of an imaginary Earth moving in the undis- 
turbed initial osculating ellipse 

equatorial radius of the central body 

difference between sines of the angle N0 at time t and 
initial time t^ for N = 1, 2, 3 

time 

initial time of osculation (epoch) 

time of periapsis passage 

argument of latitude 

argument of latitude of the Earth 

angular velocity of the spacecraft in its orbit 

angular velocity of the Earth in its orbit 

magnitude of the velocity vector of the spacecraft 

magnitude of the velocity vector of the spacecraft at time t^ 

velocity vector of the spacecraft 

velocity vector of the Earth 

potential of the central body 
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x, y, z 


x, y, z 


X E’ y E’ Z E 


X E’ y E’ Z E 


X, Y, Z 


X, Y, Z 


X E> Y E’ Z E 


X E’ Y E’ Z E 


tn= 6(°)' 


potential o£ the central body at time t^ 

heliocentric rectangular coordinates of the spacecraft in 
the equatorial reference frame of the central body 

velocity components of the spacecraft along the axes of the 
equatorial reference frame of the central body 

heliocentric rectangular coordinates of the Earth in the 
equatorial reference frame of the central body 

velocity components of the Earth along the axes of the 
equatorial reference frame of the central body 

heliocentric rectangular coordinates of the spacecraft in 
the equatorial reference frame of the Earth 

velocity components of the spacecraft along the axes of the 
equatorial reference frame of the Earth 

heliocentric rectangular coordinates of the Earth in the 
equatorial reference frame of the Earth 

velocity components of the Earth along the axes of the 
equatorial reference frame of the Earth 

right ascension of the north pole of the central body 

vernal equinox point for 1950. 0 

declination of the north pole of the central body 

a constant 

obliquity of the ecliptic 

argument of Legendre polynomials; Sine of the latitude angle 

value of £ at time t^ 

true anomaly angle of the spacecraft 

true anomaly angle of an imaginary particle moving in the 
unperturbed initial reference ellipse 

true anomaly angle of the Earth 

angular velocity of the spacecraft in its orbit 
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angular velocity of the Earth in its orbit 
true anomaly of the spacecraft at time t^ 
gravitational constant of the central body 
geocentric position vector of the observer, £ = 1 1"| 
velocity vector of the observer 
topocentric range of the spacecraft 

topocentric range rate: radial velocity of the spacecraft 

relative to the observer’s position 

topocentric position vector of the spacecraft 

velocity vector of the spacecraft relative to the observer 

topocentric position vector of an imaginary particle moving 
in the unperturbed initial ellipse 

velocity vector, relative to the observer, of an imaginary 
particle moving in the initial unperturbed orbit 

difference between the mean anomalies of the spacecraft and 
an imaginary particle moving in the undisturbed orbit 

latitude (declination) of the spacecraft 

eccentric anomaly angle of the spacecraft 

argument of perihelion of the spacecraft measured from the 
intersection of its orbital plane with the equatorial plane of 
the central body 

argument of perihelion of the Earth measured from the 
intersection of the Earth’s orbital plane with the equatorial 
plane of the central body 

argument of perihelion of the spacecraft measured from the 
intersection line of the Earth’s orbital and the equatorial 
planes 
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argument of perihelion of the Earth measured from the 
intersection line of the Earth's orbital and equatorial 
planes for 1950. 0 

argument of perihelion of the Earth measured from the 
intersection of the Earth's orbital plane and the ecliptic 
plane for 1 590 .0 

longitude of the ascending node of the spacecraft's orbital 
plane with respect to the equatorial plane of the central body- 

longitude of the ascending node of the Earth's orbital plane 
with respect to the equatorial plane of the central body 

longitude of the ascending node of the orbital plane of the 
spacecraft relative to the equatorial plane of the Earth 
for 1950.0 

longitude of the ascending node of the orbital plane of the 
Earth relative to its equatorial plane for 1950. 0 

longitude of the ascending node of the orbital plane of the 
Earth relative to the ecliptic plane for 1950. 0 
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P(t) 



Fig. 1. Initial osculating orbit 



Fig. 2. Orbital geometry 
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Fig. 3. Argument of latitude 



Fig. 4. Fundamental reference planes 
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INPUT : 

*Q' *TD/ I = ^ t • • • r n 

a(J)# e(J), i'(J) 
ft'(J), u'U) 

J = 1 SPACECRAFT 
J = 2 EARTH 


COMPUTE ANGLES 
GROUP C63 


CONSTANTS 
GROUP C2] 


i 





BB 

NO 




J = J • 


POSITION IN REFERENCE 
ELLIPSE 

r*(l,J), 9 *(1/ J)/ r*(l/J) 


Aa(l, J), Ae(l, J), 

m,j), 

Aw(U), At(I,J) 
GROUP [3] 



PARTIALS: 

DAJ(I,J), DEJ(I,J), 
DIJ(LJ), DOJ(l,J), 
DWJ(I , J), DTJ(I,J) 
GROUP \jQ 


Ar(l,J), Ar(l, J) 
Aud.j) = /-*• 
Au(l,J) = i 
GROUP GO 


Au(l,J) = Ae(l,J) + Aa>(l, J) 
Au(i / j)_=_Aed,J) 


YES j 


EQUATORIAL COORDINATES 
AND VELOCITIES 
X,Y,Z,X,Y,2(I,J) 

. GIVEN 


v NO 


TABLE OF VALUES 
X,Y # Z / ^ / Y # Zd # J) 


COMPUTE # 
x # y,z,x / y / zd,J) / 
Eqs. (73) 


COMPUTE 

x,y,z,x 

GROUP 





PARTIALS OF r,r,u,u: 
DRJ(IJ), DUJ(U) 
DDRJ(|, J), DDUJ(U) 
GROUP C123 



PARTIALS OF COORDINATES 
AND VELOCITIES 
DXJ(I,J), DYJ(I, J), DZJ(U) 
DDXJ(I, J), DDYJ(I,J), DDZ J(l, J) 
GROUPCin 


Fig. 7. Flow chart of the computational logic 
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